#include "../include/Multigrid_Solver.h"
#include "../include/Multigrid_Solver_1D.h"
#include "../include/Multigrid_Solver_2D.h"
#include "../include/Multigrid_Factory.h"
#include "../include/Functions_test.cpp"
#include <string>

std::string Rtype2str(int i)
{
    if (i == 1)
        return "FullWeight";
    else if (i == 2)
        return "Injection";
    else
        return "unkown";
};

std::string Itype2str(int i)
{
    if (i == 1)
        return "Linear";
    else if (i == 2)
        return "Quadratic";
    else
        return "unkown";
};

std::string Stype2str(int i)
{
    if (i == 1)
        return "WTJacobi";
    else if (i == 2)
        return "GaussSeidal";
    else
        return "unkown";
};

int main()
{
    int maxium_iteration = 50;
    std::cout << "Nonhomogeneous Possion Equation:" << std::endl;
    std::cout << "u''(P) =  (-2*sin(P[0])*sin(P[1])+sin(P[1])*cos(P[0])*cos(P[0])+sin(P[0])*cos(P[1])*cos(P[1]))*exp(sin(P[0])*sin(P[1])), and u=exp(sin(P[0])*sin(P[1])) on the boundary." << std::endl;
    int numofgrid;   
    double _a[2] = {0,0};
    double h;
    VectorXd u_, f_, r_, solution_temp, rhs, r_temp, r_origin, solution_exact, sol_err_box, res_err_box;
    double r_origin_norm;
    Factory1DMGSolver Solver1D;
    Factory2DMGSolver Solver2D;
    FactoryMGSolverBase *_p_fac_base = &Solver1D;
    Multigrid_Solver *_arena1D = _p_fac_base->CreateMGSolver();
    _p_fac_base = &Solver2D;
    Multigrid_Solver *_arena2D = _p_fac_base->CreateMGSolver();
    int n;
    int _Rtype;
    int _Itype;

    for (n = 128; n <= 1024; n *= 2)
    {
        for (int _Stype = 2; _Stype <= 2; _Stype++)
        {
            for (int _Rtype = 1; _Rtype <= 2; _Rtype++)
            {
                for (int _Itype = 1; _Itype <= 2; _Itype++)
                {
                    _arena2D->Initialize(n, Dim2_type1_df2rd, Dim2_type1_f);
                    _arena2D->setRestircType(_Rtype);
                    _arena2D->setIntpType(_Itype);
                    _arena2D->setSmoothType(_Stype);
                    solution_exact.setZero((n+1)*(n+1));
                    numofgrid=n;
                    for (int j = 0; j <= numofgrid; j++)
                    {
                        for (int i = 0; i <= numofgrid; i++)
                        {
                            _a[0] = 1.0 * i / numofgrid;
                            _a[1] = 1.0 * j / numofgrid;
                            solution_exact(j * (numofgrid + 1) + i) = Dim2_type1_f(_a);
                        };
                    };
                    u_ = _arena2D->get_u_origin();
                    f_ = _arena2D->get_f_origin();
                    h = _arena2D->get_h();
                    r_origin = _arena2D->ResidualSolver(u_, f_, h);
                    r_origin_norm = r_origin.lpNorm<Eigen::Infinity>();
                    u_ = _arena2D->FMGcycle(f_, h);
                    r_ = _arena2D->ResidualSolver(u_, f_, h);
                    std::cout << "FMG :n= " << n << " " <<"Type1:nonhomo" <<std::endl;
                    printf("#Smoother               #Restricter             #Interpolater \n");
                    printf("%10s    %20s    %20s\n", Stype2str(_Stype).c_str(), Rtype2str(_Rtype).c_str(), Itype2str(_Itype).c_str());
                    std::cout << "#Detail:" <<std::endl;
                    std::cout << "relative accu= " << r_.lpNorm<Eigen::Infinity>() / r_origin_norm << std::endl;
                    std::cout << "Residual norm=" << r_.lpNorm<Eigen::Infinity>() << std::endl;
                    std::cout << "Error norm=" <<(u_-solution_exact).lpNorm<Eigen::Infinity>()<<std::endl;
                    std::cout <<"Ratio(error)="<<(u_-solution_exact).lpNorm<Eigen::Infinity>()/(-solution_exact).lpNorm<Eigen::Infinity>()<<std::endl;
                    std::cout<<" "<<std::endl;
                };
            };
        };
    };
    for(int i=1;i<=6;i++)
    std::cout<<" "<<std::endl;
    std::cout << "Homogeneous Possion Equation:" << std::endl;
    std::cout << "u'' = 0, and u=x+y on the boundary." << std::endl;

        for (n = 128; n <= 1024; n *= 2)
    {
        for (int _Stype = 2; _Stype <= 2; _Stype++)
        {
            for (int _Rtype = 1; _Rtype <= 2; _Rtype++)
            {
                for (int _Itype = 1; _Itype <= 2; _Itype++)
                {
                    _arena2D->Initialize(n, Dim2_type2_df2rd, Dim2_type2_f);
                    _arena2D->setRestircType(_Rtype);
                    _arena2D->setIntpType(_Itype);
                    _arena2D->setSmoothType(_Stype);
                    solution_exact.setZero((n+1)*(n+1));
                    numofgrid=n;
                    for (int j = 0; j <= numofgrid; j++)
                    {
                        for (int i = 0; i <= numofgrid; i++)
                        {
                            _a[0] = 1.0 * i / numofgrid;
                            _a[1] = 1.0 * j / numofgrid;
                            solution_exact(j * (numofgrid + 1) + i) = Dim2_type2_f(_a);
                        };
                    };
                    u_ = _arena2D->get_u_origin();
                    f_ = _arena2D->get_f_origin();
                    h = _arena2D->get_h();
                    r_origin = _arena2D->ResidualSolver(u_, f_, h);
                    r_origin_norm = r_origin.lpNorm<Eigen::Infinity>();
                    u_ = _arena2D->FMGcycle(f_, h);
                    r_ = _arena2D->ResidualSolver(u_, f_, h);
                    std::cout << "FMG :n= " << n << " " <<"Type2:Homo"<< std::endl;
                    printf("#Smoother               #Restricter             #Interpolater \n");
                    printf("%10s    %20s    %20s\n", Stype2str(_Stype).c_str(), Rtype2str(_Rtype).c_str(), Itype2str(_Itype).c_str());
                    std::cout << "#Detail:" <<std::endl;
                    std::cout << "relative accu= " << r_.lpNorm<Eigen::Infinity>() / r_origin_norm << std::endl;
                    std::cout << "Residual norm=" << r_.lpNorm<Eigen::Infinity>() << std::endl;
                    std::cout << "Error norm=" <<(u_-solution_exact).lpNorm<Eigen::Infinity>()<<std::endl;
                    std::cout <<"Ratio(error)="<<(u_-solution_exact).lpNorm<Eigen::Infinity>()/(-solution_exact).lpNorm<Eigen::Infinity>()<<std::endl;
                    std::cout<<" "<<std::endl;
                };
            };
        };
    };
std::cout<<"EOF"<<std::endl;
};